
% Make sure to load dataFaces, dataScenes, dataBodies and dataObjects into
% workspace before executing this script


%% some constants and label definitions

% categories to be included in analysis
lab.categories            = {'Faces','Scenes','Bodies','Objects'};

% constants
con.Fs                    = 300;
con.nFacesTrials          = size(dataFaces.trial,1);
con.nScenesTrials         = size(dataScenes.trial,1);
con.nBodiesTrials         = size(dataBodies.trial,1);
con.nObjectsTrials        = size(dataObjects.trial,1);

con.timeWindowInSeconds   = 0.040;
con.timeWindowInSamples   = round(con.timeWindowInSeconds * con.Fs)-1;
con.nTimepoints           = size(dataFaces.time,2)-con.timeWindowInSamples;

% labels
lab.timeAxis              = dataFaces.time(1:con.nTimepoints);
lab.designMatrix          = [ones(con.nFacesTrials,1); 2*ones(con.nScenesTrials,1);...
        3*ones(con.nBodiesTrials,1); 4*ones(con.nObjectsTrials,1)];

    
%% cross-validation

clear decodingResults

for t = 1:con.nTimepoints
    
    display(['computing timepoint ' int2str(t) ' of ' int2str(con.nTimepoints)])
    
    squeezedDataFaces   = squeeze(mean(dataFaces.trial(:,:,t:t+con.timeWindowInSamples),3));
    squeezedDataScenes  = squeeze(mean(dataScenes.trial(:,:,t:t+con.timeWindowInSamples),3));
    squeezedDataBodies  = squeeze(mean(dataBodies.trial(:,:,t:t+con.timeWindowInSamples),3));
    squeezedDataObjects = squeeze(mean(dataObjects.trial(:,:,t:t+con.timeWindowInSamples),3));
    
    xData       = [squeezedDataFaces; squeezedDataScenes; squeezedDataBodies; squeezedDataObjects];
    
    classifier  = dml.glmnet('lambda',1e-4,'alpha',0,'family','multinomial');
    
    c           = dml.crossvalidator('mva',{dml.standardizer classifier},'stat','accuracy');
    
    c           = c.train(xData,lab.designMatrix);
    
    decodingResults.accuracy(t)  = c.statistic;
    
end


%% visualize results

plot(lab.timeAxis,decodingResults.accuracy,'r','LineWidth',2)
title('accuracy')

